function [varargout]=Equ_HH(PP,SS,MODEL,JacobianFlag,Dev_Input,Dev_Output,TimePrintFlag)

%% Preliminaries
VAR         =   MODEL.VAR;
EQ          =   MODEL.EquBlock.HH;

if nargin<=3
    JacobianFlag    =   0;
end
if nargin<=4
    Dev_Input       =   zeros(EQ.Dim.In,1);
    Dev_Output      =   zeros(EQ.Dim.Out,1);
elseif ( isempty(Dev_Input) || isempty(Dev_Output) )
    Dev_Input       =   zeros(EQ.Dim.In,1);
    Dev_Output      =   zeros(EQ.Dim.Out,1);
end
if nargin<=6
    TimePrintFlag   =   0;
end
if TimePrintFlag
    TimeStart       =   tic;
end
%% Evaluation
%--------------------------------------------------------------------------
% Pre-Transformation: Input Deviations to Input Levels
%--------------------------------------------------------------------------

for ii=1:EQ.Num.In
    TempVar     =   EQ.Var.In{ii};
    TempDev     =   Dev_Input(EQ.LocIdx.In.(TempVar));
    TempLevel   =   VAR.(TempVar).Dev2Level(TempDev);
    eval(['Dev_',TempVar,'=TempDev;'])
    eval([TempVar,'=TempLevel;'])
end
Dev_Decision=   [Dev_Pir;Dev_dE;Dev_Lag_PortfolioInfo;Dev_PortfolioInfo;...
                 Dev_T;Dev_Profit;Dev_TAU;Dev_ir_dom;Dev_ir_ext;Dev_w_N;Dev_w_H;...
                 Dev_ValFunp];
Dev_Distribution ...
            =   [Dev_Dist];

%--------------------------------------------------------------------------
% 1. Decision Problem
%--------------------------------------------------------------------------
if max(abs(Dev_Decision))<eps && JacobianFlag>0
    ValFun      =   VAR.ValFun.SS;
    Policy      =   SS.DistPolicy;
    FlowVec     =   SS.TrProb.FlowVec;
    CoefMat     =   SS.TrProb.CoefMat;
    UnitTrProb  =   SS.TrProb.UnitTrMat;
else
    %----------------------------------------------------------------------
    % Solve the Decision Problem Sequentially
    %----------------------------------------------------------------------
    % Start Step: inverse State Reduction, ValFunp to EVCoefp
    UCoefp      =   SS.UCoef+ValFunp;
    % Step 1: from UCoefp to VbarCoef
    Opt_Vbar    =   Fun_Vbar(PP,SS,T,Profit,TAU,ir_dom,ir_ext,w_N,w_H,PortfolioInfo,UCoefp,SS.FunApp.V.State);
    VbarCoef_Mat=   SS.FunApp.V.BasMat_0\reshape(Opt_Vbar,SS.FunApp.V.MatSize_State);
    VbarCoef    =   VbarCoef_Mat(:);
    % Step 2: from VCoef to VbarCoef
    Opt_U       =   Fun_U(PP,SS,Pir,dE,Lag_PortfolioInfo,VbarCoef,SS.FunApp.U.State);
    UCoef_Mat   =   SS.FunApp.U.BasMat_0\reshape(Opt_U,SS.FunApp.U.MatSize_State);
    UCoef       =   UCoef_Mat(:);
    % End Step: UCoef to ValFun
    ValFun      =   UCoef-SS.UCoef;
    %----------------------------------------------------------------------
    % Solve the Policy over Distribution Grid and Transition Prob.
    %----------------------------------------------------------------------
    % Policy over the Histogram Grid
    Policy      =   Policy_Solver(PP,SS,T,Profit,TAU,ir_dom,ir_ext,w_N,w_H, ...
                                  PortfolioInfo,UCoefp,SS.DistApp.fzo.State);
    % Transition Prob. Matrix determined by Policy
    [FlowVec,CoefMat,UnitTrProb] ...
                =   Distribution_TrProbMat(PP,SS,Pir,dE,Lag_PortfolioInfo,VbarCoef,Policy);  
    
end

%--------------------------------------------------------------------------
% 2. Distribution Evolution
%--------------------------------------------------------------------------
if max(abs(Dev_Distribution))<eps && JacobianFlag>0
    QW_U        =   SS.Dist_fzo.QW_fzo;
else
    % inverse State Reduction: (Marginal_f,Marginal_zo) to (f,zo)
    QW_U        =   DimReduction_Distribution(PP,SS,'UnZip',Dist,SS.CopulaF);
end

% Evolution: U2Up
QW_Up       =   CoefMat*QW_U+FlowVec;
% State Reduction:
[Distp,~]   =   DimReduction_Distribution(PP,SS,'Zip',QW_Up);

%--------------------------------------------------------------------------
% 3. Aggregation
%--------------------------------------------------------------------------
Agg         =   Aggregation(PP,SS,QW_U,Policy,UnitTrProb);

Bp_dom      =   Agg.Bp_dom;

B_dom       =   Agg.B_dom;
B_ext       =   Agg.B_ext;

C           =   Agg.C;
C_FI_dom    =   Agg.Tot.c.dom;
C_FI_ext    =   Agg.Tot.c.ext;
C_RI_N      =   Agg.Tot.c.N;
C_RI_H      =   Agg.Tot.c.H;

L_N         =   Agg.L_N;
L_H         =   Agg.L_H;

AggStat     =   [Agg.Avg.c.dom;Agg.Avg.c.ext;Agg.Avg.c.N;Agg.Avg.c.H;Agg.Avg.c.poor;Agg.Avg.c.rich;...
                 Agg.Avg.c.dom_N;Agg.Avg.c.dom_H;Agg.Avg.c.ext_N;Agg.Avg.c.ext_H;...
                 Agg.Avg.c.dom_poor;Agg.Avg.c.dom_rich;Agg.Avg.c.ext_poor;Agg.Avg.c.ext_rich;...
                 Agg.Avg.c.H_poor;Agg.Avg.c.H_rich;Agg.Avg.c.N_poor;Agg.Avg.c.N_rich];
%--------------------------------------------------------------------------
% Post-Transformation: Output Levels to Ouput Deviations
%--------------------------------------------------------------------------
Dev_Left    =   zeros(EQ.Dim.Out,1);
for ii=1:EQ.Num.Out
    TempVar     =   EQ.Var.Out{ii};
    eval(['TempLevel=',TempVar,';']);
    TempDev     =   VAR.(TempVar).Level2Dev(TempLevel);
    Dev_Left(EQ.LocIdx.Out.(TempVar)) ...
                =   TempDev;
end
%% Return Objectives

switch JacobianFlag
    case 0
        varargout{1}    =   -Dev_Output+Dev_Left;
    case 1
        TempFun         =   @(Dev_Input)Equ_HH(PP,SS,MODEL,0,Dev_Input,Dev_Output);
        Jac             =   struct();
        Jac.Input       =   parfdjac(TempFun,Dev_Input);
        Jac.Output      =   -speye(MODEL.EquBlock.HH.Dim.Out);
        
        for ii=1:EQ.Num.In
            TempVar         =   EQ.Var.In{ii};
            Jac.(TempVar)   =   Jac.Input(:,EQ.LocIdx.In.(TempVar));
        end
        for ii=1:EQ.Num.Out
            TempVar         =   EQ.Var.Out{ii};
            Jac.(TempVar)   =   Jac.Output(:,EQ.LocIdx.Out.(TempVar));
        end
        varargout{1}    =   Jac;
end
if TimePrintFlag
    TimeElapsed     =   toc(TimeStart);
    fprintf(['Jacobian of HH Block is done: ',...
             'Elapse Time=',num2str(TimeElapsed,'%.2f'),'s','\n']);
end

